### Account for uncertainty in latent variable measure ###
##########################################################

# Create dataset. Users should normally import their own data.
dat <- data.frame(
c(-1.38098262083418,-1.45391602206088,-1.45273021607724,-1.44056715623687,-1.47659229832128),
c(0.30054733394201,0.25383554283216,0.25416831765569,0.27464937781364,0.2659630134477),
c(-1.45391602206088,-1.45273021607724,-1.44056715623687,-1.47659229832128,-1.50500598714741),
c(0.25383554283216,0.25416831765569,0.27464937781364,0.2659630134477,0.26203514598166)
)
colnames(dat) <- c("latentmean","latentsd","latentmean.lag","latentsd.lag")

# Read in the dataset
# data <- read.csv()

# Create a list to save datasets
newdata <- list()

# Create a data frame that includes your model variables
datnew <- na.omit(data.frame(dat$latentmean,dat$latentsd,dat$latentmean.lag,dat$latentsd.lag))

# Take K draws from the posterior distribution and make K new datasets in the list object just created
for(i in 1:10){
  newdata[[i]] <- datnew
  newdata[[i]]$latentmean <- rnorm(cbind(rep(1,nrow(dat))), mean=datnew$dat.latentmean, sd=datnew$dat.latentsd)
  newdata[[i]]$latentmean.lag <- rnorm(cbind(rep(1,nrow(dat))), mean=datnew$dat.latentmean.lag, sd=datnew$dat.latentsd.lag)
  
}

# Inspect list object. NOTE: Excercise caution when inspecting the list object while using a full dataset and/or for large k.
newdata

# Define a formula as a character string like this: FML <- "Y ~ X"
FML <- "latentmean ~ latentmean.lag"

# Create function. The code below runs a linear model K times and combines all the estimates using the formula for combining multiply imputed datasets from Rubin.
milm <- function(fml, midata){
  xx <- terms(as.formula(fml))
  lms <- matrix(data=NA, nrow=(length(attr(xx, "term.labels")) + 1), ncol=length(midata))
  ses <- matrix(data=NA, nrow=(length(attr(xx, "term.labels")) + 1), ncol=length(midata))
  vcovs <- list()
  for(i in 1:length(midata)){
    tmp <- lm(formula=as.formula(fml), data=midata[[i]])
    lms[,i] <- tmp$coefficients
    ses[,i] <- sqrt(diag(vcov(tmp)))
    vcovs[[i]] <- vcov(tmp)
  }
  par.est <- apply(lms, 1, mean)
  se.within <- apply(ses, 1, mean)
  se.between <- apply(lms, 1, var)
  se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(midata))))
  list("terms"=names(tmp$coefficients), "beta" = par.est, "SE"=se.est, "vcovs"=vcovs,"coefs" = lms)   
}

# Call function
model <- milm(fml=FML, midata=newdata)

# Inspect returned object
attributes(model)

# Organize results
tstats <- data.frame(model$terms,model$beta,model$SE,model$beta/model$SE)
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat")

# View results
tstats